﻿using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using System.Threading.Tasks;

namespace CulvertDesigner
{

    public static class Bisect
    {
        public static double FindRoot(Func<double, double> func, double a, double b, double tol = 1e-5, int maxIter = 100)
        {
            double fa = func(a);
            double fb = func(b);

            if (Math.Sign(fa) == Math.Sign(fb))
            {
                throw new ArgumentException("The function values at the interval endpoints must have opposite signs.");
            }

            double c = a;
            for (int i = 0; i < maxIter; i++)
            {
                double fc = func(c);
                if (Math.Abs(fc) < tol)
                {
                    return c;
                }

                c = (a + b) / 2;
                double fc_new = func(c);

                if (Math.Sign(fc_new) == Math.Sign(fa))
                {
                    a = c;
                    fa = fc_new;
                }
                else
                {
                    b = c;
                    fb = fc_new;
                }
            }

            throw new InvalidOperationException("Bisection method failed to converge.");
        }
    }

    public static class CulvertCalculator
    {
        /// <summary>
        /// 计算圆截面流体面积
        /// </summary>
        /// <param name="diameter">直径</param>
        /// <param name="theta">角度</param>
        /// <returns>圆截面流体面积</returns>
        public static double CalculateArea(double diameter, double theta)
        {
            // 使用公式计算圆截面流体面积
            return (theta - Math.Sin(theta)) / 8 * Math.Pow(diameter, 2);
        }

        /// <summary>
        /// 根据深度计算圆心角
        /// </summary>
        /// <param name="depth">深度</param>
        /// <param name="diameter">直径</param>
        /// <returns>圆心角</returns>
        public static double CalculateThetaByDepth(double depth, double diameter)
        {
            double rr = diameter * 0.5;
            if (depth / diameter > 1)
            {
                throw new ArgumentException("Depth cannot exceed diameter.");
            }
            else
            {
                double theta;
                if (depth < rr)
                {
                    theta = 2.0 * Math.Acos((rr - depth) / rr);
                }
                else
                {
                    theta = 2.0 * Math.PI - 2.0 * Math.Acos((depth - rr) / rr);
                }
                return theta;
            }
        }

        /// <summary>
        /// 计算圆截面流体深度
        /// </summary>
        /// <param name="diameter">直径</param>
        /// <param name="theta">角度</param>
        /// <returns>圆截面流体深度</returns>
        public static double CalculateDepth(double diameter, double theta)
        {
            // 使用公式计算圆截面流体深度
            return 0.5 * (1 - Math.Cos(theta / 2)) * diameter;
        }


        /// <summary>
        /// 计算圆截面流体湿周长
        /// </summary>
        /// <param name="diameter">直径</param>
        /// <param name="theta">角度</param>
        /// <returns>圆截面流体湿周长</returns>
        public static double CalculatePerimeter(double diameter, double theta)
        {
            // 使用公式计算圆截面流体湿周长
            return theta / 2 * diameter;
        }

        /// <summary>
        /// 计算圆截面流体宽度
        /// </summary>
        /// <param name="diameter">直径</param>
        /// <param name="theta">角度</param>
        /// <returns>圆截面流体宽度</returns>
        public static double CalculateWidth(double diameter, double theta)
        {
            // 使用公式计算圆截面流体宽度
            return Math.Sin(theta / 2) * diameter;
        }

        public static double CalculateError1(double th, double dia, double discharge)
        {
            return Math.Pow(discharge, 2) / 9.806 - Math.Pow(CalculateArea(dia, th), 3) / CalculateWidth(dia, th);
        }

        public static double CalculateError2(double th, double dia, double flow, double s0, double n)
        {
            double wp = CalculatePerimeter(dia, th);
            double area = CalculateArea(dia, th);
            double rh = area / wp;
            double v = 1.0 / n * Math.Pow(rh, 2.0 / 3) * Math.Pow(s0, 0.5);
            return area * v - flow;
        }

        public static double CalculateUniformDepth(double flow, double diameter, double s0, double n)
        {
            double wp = Math.PI * diameter;
            double area = Math.Pow(diameter, 2) * 0.25 * Math.PI;
            double rh = area / wp;
            double vmax = 1.0 / n * Math.Pow(rh, 2.0 / 3) * Math.Pow(s0, 0.5);
            double Qmax = vmax * area;

            if (flow >= Qmax)
            {
                return diameter;
            }
            else
            {
                Func<double, double> errorFunction = th => CalculateError2(th, diameter, flow, s0, n);
                try
                {
                    double r = Bisect.FindRoot(errorFunction, 1e-8 * Math.PI, 2 * Math.PI);
                    return CalculateDepth(diameter, r);
                }
                catch (Exception ex)
                {
                    throw new InvalidOperationException("Error occurred during depth calculation.", ex);
                }
            }
        }
        public static double CalculateErrorH(double qTest, double hTar, double dia, double s0, double n)
        {
            return CalculateUniformDepth(qTest, dia, s0, n) - hTar;
        }

        public static (double, string) CalculateCapacity(double diameter, double s0, double n, double limH = 0.75, double limV = 5.0)
        {
            Func<double, double> errorFunction = qTest => CalculateErrorH(qTest, limH * diameter, diameter, s0, n);
            double maxQByH = Bisect.FindRoot(errorFunction, 1e-3, 12);

            // 如果需要计算 max_q_by_v，取消下面的注释并实现 _error_v 函数的转换
            // Func<double, double> errorVFunction = vTest => ErrorVCalculator.CalculateError(vTest, limV, diameter, s0, n);
            // double maxQByV = Bisect.FindRoot(errorVFunction, 1e-3, 12);

            return (maxQByH, $"{limH:F2}D");
        }

        /// <summary>
        /// 计算稳流的福禄数
        /// </summary>
        /// <param name="flow">流量</param>
        /// <param name="diameter">管径</param>
        /// <param name="s0">底坡</param>
        /// <param name="n">Manning系数</param>
        /// <returns>福禄数</returns>
        public static double CalculateFroudeNumber(double flow, double diameter, double s0, double n)
        {
            // 计算均匀流深度
            double du = CalculateUniformDepth(flow, diameter, s0, n);

            // 根据深度计算角度
            double th = CalculateThetaByDepth(du, diameter);

            // 根据角度计算宽度
            double T = CalculateWidth(diameter, th);

            // 根据深度计算面积
            double A = CalculateArea(diameter, th);

            // 根据角度计算湿周长
            double wp = CalculatePerimeter(diameter, th);

            // 计算平均水深
            double dm = A / T;

            // 计算流速
            double v = flow / A;

            // 计算福禄数
            return v / Math.Sqrt(9.806 * dm);
        }

        /// <summary>
        /// 计算圆截面的临界深度
        /// </summary>
        /// <param name="Q">流量</param>
        /// <param name="diameter">直径</param>
        /// <returns>临界深度</returns>
        public static double CalculateCriticalDepth(double Q, double diameter)
        {
            Func<double, double> errorFunction = th => CalculateError1(th, diameter, Q);
            try
            {
                double r = Bisect.FindRoot(errorFunction, 1e-3 * Math.PI, 2 * Math.PI);
                return CalculateDepth(diameter, r);
            }
            catch (Exception ex)
            {
                throw new InvalidOperationException("Error occurred during critical depth calculation.", ex);
            }
        }

        /// <summary>
        /// 计算临界流速
        /// </summary>
        /// <param name="Q">流量</param>
        /// <param name="diameter">直径</param>
        /// <returns>临界流速</returns>
        public static double CalculateCriticalVelocity(double Q, double diameter)
        {
            double dc = CalculateCriticalDepth(Q, diameter);
            double a = CalculateArea(diameter, CalculateThetaByDepth(dc, diameter));
            return Q / a;
        }


    }




}
